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1 .  Introduction 


The  research  on  Stochastic  and  Centrifuge  Modelling  of  Jointed 
Rock  has  been  performed  during  the  years  1987  to  1990  under 
sponsorship  of  the  U.S.  Air  Force  Office  of  Scientific  Research  and  the 
U.S.  Air  Force  Engineering  Services  Center.  It  represents  a  logical 
continuation  of  preceding  work  by  the  Pi's  on  geometric  and 
mechanical  modelling  of  jointed  rock  and  centrifuge  studies  on 
granular  material. 

The  research  was  subdivided  into  three  major  areas,  namely: 

Fracturing  of  Jointed*  Rock 

Centrifuge  Modelling  of  Jointed  Rock 

Stochastic  Modelling  of  Jointed  Rock 

which  are  parts  of  an  entity,  the  behavior  of  a  jointed  mass. 

Stochastic  modelling  of  jointed  rock  is  necessary  to  represent  the 
geometry  of  rock  masses.  Fracturing  of  jointed  rock  is  a  possible 
mechanism  underlying  rock  mass  behavior  in  which  a  jointed  (or 
discontinuous  mass  deforms  and  eventually  fails  as  new  fractures 
connect  existing  ones.  The  other  major  mechanism  of  rock  mass 
failure  is  relative  sliding  of  blocks  bounded  by  existing  joints  or  by 
existing  joints  and  newly  created  fractures;  this  can  be  physically 
modelled  in  the  centrifuge  under  realistic  stress  conditions. 

The  combination  of  geometric  and  mechanical  modelling  of  a 
jointed  mass  seems  to  be  an  obvious  approach.  Previous  work  of  our 
group  and  others  in  the  area  suffer,  however,  from  several 
inadequacies  or  lack  of  knowledge.  Stochastic  and  other  geometric 
models  which  existed  prior  to  the  research  reported  here  were 
usually  not  capable  of  representing  the  actual  geometry  of  a  rock 
mass  even  if  the  information  mass  was  carefully  sampled  from  it. 


*  Throughout  this  report  the  terms  "joint"  and  "fracture"  will  be  used 
interchangeably. 


2 


The  occurrence  of  new  fractures  in  a  discontinuum  was  not  well 
understood,  in  particular,  only  limited  experimental  observations  and 
analytical  models  existed.  While  the  basic  understanding  of 
movement  of  rigid  blocks  bounded  by  discontinuities  exists, 
analytical  modelling  usually  relies  on  limit  equilibrium  approaches 
which  neglect  the  real  stress  distribution.  Physical  experiments  were 
not  adequate  to  represent  geostatic  stress  distribution,  or  if  they 
were,  such  as  centrifuge  experiments,  doubts  existed  as  to  the 
applicability  of  centrifuge  testing  to  discontinuua. 

The  research  described  in  this  report  made  significant 
contributions  in  the  three  problem  areas,  through  the  development 
of  a  stochastic  model  which  can  represent  real  joint  geometry, 
through  experimentation  and  analytical  modelling  of  fracturing 
between  existing  joints,  and  through  centrifuge  experimentation 
which  establishes  the  limitations  of  this  method  with  regard  to 
testing  of  discontinua.  The  research  went  beyond  these  contributions 
in  that  the  stochastic  model  was  combined  with  a  topological  model 
for  slope  stability,  in  that  the  fracture  coalescence  work  provided 
general  insight  into  fracturing  of  discontinuous  brittle  materials,  and 
in  that  the  centrifuge  investigation  led  to  a  complete  understanding 
and  analytical  modelling  for  arching  of  granular  continua  and  of 
discontinuua.  As  we  will  show  in  the  conclusions,  this  research  not 
only  led  to  major  advances  in  the  fundamental  understanding  of 
geometric  and  mechanical  behavior  of  rock  masses,  but  it  also 
advanced  fundamental  knowledge  on  brittle  discontinua  in  general, 
and  it  provided  results  which  have  substantial  practical  significance. 

The  report  is  subdivided  into  three  parts  corresponding  to  the 
three  research  areas: 

Part  I:  Fracturing  of  Fractured  Rock 

Part  II:  Centrifuge  Modelling  of  Jointed  Rock 

Part  III:  Stochastic  and  Topological  Fracture  Geometry  Model 
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each  of  which  is  contained  in  a  separate  volume.  In  the  following 
lines  we  will  summarize  in  the  same  sequence  the  major  aspects  of 
the  work  and  the  most  important  results. 
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2.  Fracturing  of  Fractured  Rock 

2.1  Problem  Statements  and  Schedule  of  Research 

Fracturing  of  'intact'  (or  otherwise  expressed,  continuous) 
brittle  material  has  received  attention  over  many  years  both  in 
theoretical  and  experimental  research.  Much  less  work  has  been 
performed  on  fracturing  of  fractured  rock,  i.e.  of  the  rock  mass 
containing  discontinuities.  This  is  so  inspite  of  the  fact  that  many 
natural  and  engineering  problems  associated  with  rock  masses 
involve  the  propagation  of  fractures  between  existing  discontinuities 
such  as  faults,  joints,  foliation  surfaces  and  bedding  planes.  Examples 
where  this  process  is  important,  are  natural  phenomena  such  as  joint 
genesis,  faulting  and  landslides  and  engineered  structures  such  as 
slopes  and  tunnels.  Although  usually  on  a  smaller  scale,  similar 
fracturing  between  existing  fractures  or  cracks  occurs  in  other 
materials,  notably,  concrete  but  also  in  polymers  and  ice.  What  is 
discussed  here  under  the  catchphrase  "fracturing  of  fractured  rock" 
includes  thus  all  problems  where  a  'new'  fracture  emanates  from  an 
existing  discontinuity  and  propagates  to  connect  with  another 
discontinuity. 

During  the  first  year  of  research  (see  Annual  Report  1987/88) 
the  research  was  divided  into  two  major  areas,  namely,  the 
development  of  a  numerical  model  for  tension  softening  materials 
and  the  preparation  and  performance  of  initial  tests  on  fracturing 
model  materials  with  pre-existing  fractures.  The  numerical  model 
for  a  tension  softening  material  was  not  only  developed  but  also 
validated.  All  this  is  reported  in  addition  to  the  Annual  Report 
87/88,  in  the  M.S.  thesis,  "Numerical  Modelling  of  Fracture 
Propagation  in  Tension  Softening  Materials"*  ,  by  O.M.  Reyes  and  a 
paper,  "Numerical  Modelling  of  Fracture  Propagation  in  Tension 
Softening  Materials",  by  Reyes,  Li,  Einstein  which  is  ready  to  be 
submitted  to  a  professional  journal.  In  the  experimental  work 
during  the  first  year,  we  selected  the  model  materials,  a  gypsum  and 
an  epoxy  resin.  The  material  properties  were  determined,  and,  most 


*For  literature  references,  see  Appendix  A. 
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importantly,  a  variety  of  procedures  for  molding  of  pre-fractured 
model  specimens  were  investigated.  Initial  tests  were  run  which 
revealed  different  types  of  cracks  emanating  from  the  pre-existing 
fractures  or  flaws. 

These  initial  tests  and  follow-up  tests  during  the  second  year 
(see  Annual  Report  1988/89)  indicated  the  importance  of  procedures 
which  would  allow  us  to  observe  the  crack  propagation  in  detail. 
Therefore,  a  substantial  effort  was  developed  to  the  study  of 
different  observation  methods  including  high  speed  film  and 
videotaping.  The  required  resolutions  led  us  to  select  a  procedure  in 
which  observation  through  a  microscope  (~20  x  magnification)  is 
combined  with  the  possibility  to  taking  photographs.  Test  series 
conducted  with  this  equipment  then  revealed  the  surprising  fact 
that,  in  contrast  to  earlier  analytical  studies  and  experiments, 
coalescence  of  pre-existing  non-overlapping  fractures  in  gypsum 
occurred  through  so  called  secondary  cracks  which  develop  after  the 
well  known  wing  crack  (Fig.  2.1).  Simultaneously,  we  investigated 
possible  analytical  and  numerical  approaches  to  model  this  behavior 
which  led  us  to  the  conclusion  that  we  had  to  develop  our  own 
model. 

The  detailed  test  series  in  which  we  observed  fracture 
coalescence  and  the  development  of  the  associated  predictive  model 
were  the  subject  of  the  third  year  of  research  and  are  described  in 
Volume  I  of  this  report  series  as  well  as  in  the  following  summary. 

2.2  Coalescence  Experiments 

2.2.1  Procedures 


Uniaxial  compression  tests  were  run  on  prismatic  specimens  as 
shown  in  Fig.  2.2.  The  material  used  is  a  hydrated  gypsum  made 
from  a  mixture  of  gypsum  (Kydrocal  B-ll),  celite  (diaiomaceous 
earth)  and  water  at  ratios  watengypsum  =  0.4  and  water:  celite  =  35. 
This  gypsum  is  essentially  the  same  material  as  that  used  in 
preceding  MIT  model  rock  experiments  because  it  is  well  suited  for 


Observed  Crack  Types:  E.W.C.=External  Wing  Crack; 
I.W.C. internal  Wing  Crack;  S.C.=Secondary  Coalescing 
Crack  (Called  "Secondary  Crack"  in  the  text); 
S.W.C.=Secondary  Wing  Crack. 
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The  geometry  of  the  pre-existing  fractures  was  chosen  as 
follows: 

•  Each  block  had.  two  12.7  mm  long  pre-existing  fractures  separated 
by  a  rock  bridge  which  was  also  12.7  mm  in  length.  The  fractures 
had  a  finite  open  width  of  ~0.25  mm  and  stayed  open  during 
compression  testing.  This  had  the  advantage  of  eliminating 
friction  between  fracture  surfaces  which  is  difficult  to  quantify. 

•  The  fractures  were  inclined  at  30°,  45°  and  60°  to  the  horizontal 
with  both  fractures  always  at  the  same  inclination  in  a  particular 
specimen. 

•  The  rock  bridge  inclination  was  the  other  main  variable  in  this 
study  (see  Figure  2.3).  The  fracture  geometries  can  be 
characterized  as  either  overlapping  or  non-overlapping, 
depending  on  the  inclination  of  the  rock  bridge  with  respect  to  the 
direction  of  the  applied  uniaxial  compressive  load  (see  Figure  2.4) 

As  mentioned  before,  we  used  a  microscope  (20x)  to  monitor 
crack  growth.  Since  the  field  of  view  through  the  microscope  is 
rather  small,  step  loading  was  also  adopted  in  order  to  allow  enough 
time  for  the  entire  surface  to  be  scanned  after  each  load  increment. 
This  was  possible  because  crack  propagation  in  compression  is  a 
stable  process  prior  to  strength  failure,  and  no  creep  could  take  place 
during  the  short  time  intervals  that  the  load  was  held  constant  at 
each  step.  Steps  were  displacement  controlled  with  displacements  of 
approximately  0.005  to  0.015  mm  per  step. 

2.2.2  Observations  and  Results  of  Coalescence  Experiments 

For  all  fracture  geometries,  wing  cracks  at  the  outer  and  inner 
fracture  tips  appeared  first.  This  was  followed  by  fracture 
coalescence  through  wing  crack  growth  for  overlapping  fractures, 
and  coalescence  through  secondary  crack  growth  for  non-overlapping 
fractures.  Fig.  2.5  shows  coalescence  between  overlapping  30°-135° 
(30°  fracture  inclination-1350  rock  bridge  inclination)  fractures;  in 
this  case,  the  lower  internal  wing  crack  propagated  towards  the  tip  of 
the  upper  fracture.  (The  term  "internal"  and  "external"  refer  to 


f "^v 

\  rock  bridge  inclination  angle 


yy  fracture  inclination 
O'- - i- angle 


Figure  2.3.  Fracture  and  Rock  Bridge  Geometries 
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Figure  2.4.  Non-overlapping  (left)  and  Overlapping  (right)  Fractures 


Figure  2.5.  Fracture  Coalescence  for  Overlapping  Pre-existing 
Fractures 
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location  of  the  wing  crack  with  respect  to  the  pair  of  existing 
fractures,  see  Fig.  2.1)  Figs.  2.6-2.9  show  the  load-displacement 
curve  (Fig.  2.6)  and  photographs  (Fig.  2.7-2.9)  of  the  fracture 
coalescence  process  for  the  non-overlapping  30°-45°  fractures  (30° 
fracture  inclination  -  45°  rock  bridge  inclination).  For  this  geometry, 
wing  cracks  first  appeared  at  an  applied  load  of  25  kN  (equivalent  to 
11  MPa  assuming  uniform  distribution  on  the  loading  surface).  Fig. 
2.7,  a  photograph  taken  at  a  load  of  47  kN  (21  MPa),  shows  these 
wing  c-  as  well  as  some  spalling  near  the  fracture  tips.  At  an 

applied  load  of  54  kN  (24  MPa),  coalescence  suddenly  occurred  (Fig. 
2.8);  this  happened  so  rapidly  that  we  were  not  able  to  establish  if 
the  propagation  direction  was  from  the  rock  bridge  center  toward 
the  fracture  tips  or  vice  versa.  Fig.  2.9  is  an  expanded  view  of  the 
30°-45°  fractures  after  coalescence.  Note  that  the  internal  wing 
cracks  are  closed  while  the  external  ones  are  open,  perhaps 
indicating  that  the  stresses  which  were  keeping  the  internal  wing 
cracks  open  were  relieved  by  the  formation  of  the  coalescence  crack. 

The  coalescence  process  observed  with  30°-45°  pre-existing 
fractures  is  essentially  the  same  for  all  other  non-overlapping 
fractures,  i.e.,  a  secondary  crack  grows  within  the  rock  bridge  and 
connects  the  inner  fracture  tips.  Differences  exist  with  regard  to:  (1) 
measured  loads  at  which  coalescence  occurred  and  (2)  surface 
characteristics  of  coalescence  cracks  as  will  be  discussed  below. 

In  addition  to  observing  crack  geometry  and  sequence,  we  also 
studied  the  details  of  the  crack  surfaces.  This  fractographic  study  is 
described  in  detail  in  Volume  I  of  this  report  series.  It  showed  us 
that  coalescence  cracks  may  involve  a  complex  interaction  of  tensile 
and  shear  fractures  as  well  as  post  fracture  shearing. 

In  general,  for  the  same  fracture  inclination  angles,  the 
coalescence  loads  were  lower  for  non-overlapping  compared  to 
overlapping  fracture  geometries.  Non-overlapping  fractures  are  thus 
more  critical  at  least  in  unconfined  compression.  Furthermore, 
secondary  crack  growth  associated  with  the  coalescence  of  non¬ 
overlapping  fractures  cannot  be  explained  by  existing  fracture 
mechanics  theories.  We,  therefore,  concentrated  on  non-overlapping 
fractures  in  our  predictive  model  as  will  be  discussed  below  and  in 
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Figure  2.6.  Load-Displacement  Curve  for 
on  30°-4i5o  specimen  (see  also 


I’niaxial  Compression  Test 
Figs.  2.6  and  2.7) 


Level  =  47  kN,  Wing  Cracks  and  Some  Spalling  Can  Be 
Observed 


Figure  2.9.  Uniaxial  Compression  Test  on  30°-45°  Specimen,  (Same 
test  as  shown  in  Figs.  2.5-2.7),  Overall  View 
Note:  Open  external  wing  cracks,  closed  internal  wing 
cracks  and  secondary  coalescence  crack. 
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Volume  I.  Figure  2.10  shows  a  plot  of  the  observed  coalescence  loads 
versus  rock  bridge  orientation  angle  fot  non-overlapping  fractures. 
Triangular,  square  and  circular  symbols  refer  to  fracture  geometries 
with  fractures  that  are  inclined  at  30°,  45°  and  60°  with  respect  to 
the  horizontal.  The  following  can  be  deduced  from  Figure  2.10:  (1) 
measured  coalescence  loads  for  specimens  with  the  same  fracture 
geometries  are  practically  identical  (two  were  tested  for  each 
geometry)  confirming  the  repeatability  of  the  observed  coalescence 
phenomenon,  (2)  for  the  same  rock  bridge  orientation  angle, 
coalescence  loads  increased  with  increasing  fracture  inclination  angle, 
and  (3)  for  the  same  fracture  inclination  angle,  coalescence  loads 
were  lowest  for  the  rock  bridge  orientation  angle  of  75°.  These 
trends  had  to  be  reproduced  with  our  analytical  model  which  will 
now  be  discussed. 


2.3.  Analytical  Modelling  of  Fracture  Coalescence 
2.3.1  Basic  considerations  and  initial  modelling 

The  stress/strain  fields  within  the  rock  bridge  after  wing  crack 
growth  and  just  prior  to  secondary  crack  growth  had  to  be 
established  in  order  to  understand  the  mechanisms  underlying  the 
fracture  coalescence  observed  in  the  experiments.  Stress  analyses 
were  thus  performed  using  the  finite  element  method  and  assuming 
linear  elastic  conditions  (E  =  6200  MPa,  v  =  0.28).  The  latter  was  a 
reasonable  assumption  since  the  load-displacement  curves  showed 
linear  elastic  behavior  up  to  the  point  of  brittle  failure. 

A  major  finding  in  these  linear  elastic  stress  analyses  was  that 
tensile  stresses  within  the  rock  bridge  were  not  relieved  by  wing 
crack  growth,  as  can  be  seen  by  comparing  Fig.  2.11  (a)  and  (b) 
which  are  vector  plots  of  tensile  major  principal  stresses  around  30°- 
75°  fractures  before  and  after  wing  crack  growth,  respectively.  The 
tensile  stresses  in  the  rock  bridge  zone  in  Fig.  2.11  (b)  indicate  that 
crack  initiation  within  the  rock  bridge  would  not  be  surprising  at  all 
for  materials  such  as  rock  which  are  weak  in  tension. 
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Figure  2.10. 


Uniaxial  Tests  on  a  Variety  of  Specimens  with  Two 
Non-overlapping  Fractures 


Figure  2.11a,  Maximurr  Tensile  Principal  Stress  Vectors  of  Specimen 
with  30°-75°  Pre-existing  Fractures  Subject  to  Uniaxial 
Compressive  Stress  in  Vertical  Direction  (before  wing 
crack  growth). 
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Figure  2.11b.  Maximum  Tensile  Principal  Stress  Vectors  of  Specimen 
with  30°-75°  Pre-existing  Fractures  and  Internal 
Cracks  Subject  to  Uniaxial  Compressive  Stress  in 
Vertical  Direction 

Note:  High  tensile  stress  level  in  zone  between  wing 
cracks. 
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Given  the  existence  of  tensile  stresses  within  the  rock  bridges 
where  secondary  crack  growth  occurred,  it  seemed  oKvious  to  apply 
a  crack  propagation  criterion  based  on  tensile  stresses  such  as  those 
discussed  in  the  introduction.  However,  a  more  careful  study  of  the 
contour  plot  of  major  principal  stresses  around  the  30°-45°  fractures, 
shown  in  Fig.  2.12,  made  us  believe  that  a  tensile  stress-based 
criterion  was  not  appropriate  for  simulating  the  fracture  behavior  in 
the  experiments.  Specifically,  there  were  regions  of  high  tensile 
stresses  in  areas  outside  of  the  rock  bridge  where  no  cracking  was 
observed.  On  the  other  hand,  a  contour  plot  of  major  principal 
strains  (Fig.  2.13)  shows  regions  cf  high  tensile  strains  only  in  areas 
where  cracking  occurred,  such  as  close  to  the  fracture  tips  and  within 
the  rock  bridge  (see  Fig.  2.8).  Given  that  regions  of  high  tensile 
strains  better  correspond  to  the  observed  secondary  cracking,  we 
chose  to  use  the  maximum  tensile  principal  strain  rather  than  the 
maximum  tensile  principal  stress  as  a  basis  for  a  fracture  criterion. 

(It  must  emphasized  that  in  multiaxial  stress  fields  the  maximum 
principal  strain  is  not  simply  proportional  to  the  maximum  principal 
stress  in  the  same  direction.  Tensile  principal  strains  can  thus  be 
higher  for  regions  where  there  is  a  combination  of  large  tensile  and 
compressive  principal  stresses.)  It  is  interesting  to  note  that  strain- 
based  failure  criteria  can  account  for  tensile  failure  (due  to  tensile 
strains)  under  biaxial  (triaxial)  compression  even  in  a  continuum. 

This  is  in  contrast  to  tensile  failure  based  on  stress  criteria  which 
requires  the  existence  of  microcracks. 

2.3.2  Simulating  Crack  Coalescence  Using  a  Strain-based  Smeared 
Crack  Model 


In  modelling  fracture  coalescence,  two  possible  approaches 
may  be  taken:  (1)  the  discrete  crack  approach  where  the  stress  field 
around  a  crack  tip  coupled  with  a  propagation  criterion  determines 
whether  a  crack  will  lengthen  and  in  what  direction  it  would  extend 
and  (2)  the  smeared  crack  approach  where  fractures  are  assumed  to 
occur  in  locations  of  high  stress  or  strain  concentration.  We 
combined  the  smeared  crack  approach  with  a  stiffness  degradation 
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Figure  2.12. 


Maximum  Tensile  Principal  Stress  Contours  for 
Specimen  with  30°-45°  Pre-existing  Fractures  at  an 
Applied  Uniaxial  Compressive  Load  of  24  MPa.  Dotted 
Areas  Indicate  Regions  of  High  Tensile  Stress. 
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1 

Maximum  Tensile  Principal  Strain  Contours  for 
Specimen  with  30°-45°  Pre-existing  Fractures  at  an 
Applied  Uniaxial  Compressive  Load  of  24  MPa  Dotted 
Areas  Indicate  Regions  of  High  Tensile  Strain. 
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Figure  2.13. 
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occur  in  locations  of  high  stress  or  strain  concentration.  We 
combined  the  smeared  crack  approach  with  a  stiffness  degradation 
model  ("damage"  model)  .[Rots  et  al.,  1985;  Lemaitre,  1986].  This  is 
admittedly  a  simplification  of  the  actual  micromechanical  process 
which  seems  to  be  justified  at  the  present  stage  where  the 
understanding  of  the  actual  propagation  process  is  limited  (see 
Section  2.2.2). 

The  stiffness  degradation  model  which  we  have  used  to  simulate 
fracture  coalescence  is  defined  as  follows  [Simo  and  Ju,  1987]: 


^ij  =  —  D)  Qjkl£kl  (2.1) 

where  aij,  ekl  and  Qjkl  are  components  of  the  stress,  strain  and 
elastic  constitutive  tensors,  respectively,  and  D  is  an  internal  state 
variable  which  measure  the  damage  in  the  material.  D  ranges  from  0 
to  1  corresponding,  respectively,  to  a  material  having  the  original 
elastic  stiffness  (i.e.,  "no  damage")  to  one  having  zero  stiffness  (i.e., 
"completely  damaged").  In  our  approach  in  which  we  use  a  strain 
dependent  criterion  D  is  dependent  on  an  equivalent  strain  eeq, 
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where 
H(e)  =  1-exp 

and  the  equivalent  strain  is  taken  as  [Lemaitre,  1986]: 
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rate  at  which  D  increases  with  strain.  Eq.  (2.3)  also  implies  that  the 
damage  has  a  maximum  value  of  H(e0/A)  which  will  be  referred  to  as 
Dsat  (saturation  value)  in  the  subsequent  discussion. 

The  model  described  above  was  implemented  as  a  material 
subroutine  for  ABAQUS  [ABAQUS  User's  Manual,  1990].  Simulations 
were  performed  on  "edge"  models  (see  Fig.  2.14)  which  were  thought 
to  be  adequate  for  investigating  coalescence  crack  growth  in  the  rock 
bridge.  Using  edge  models  required  less  computation  time  due  to  the 
smaller  number  of  finite  elements  when  compared  to  the  "full" 
models.  (The  elastic  strain  fields  within  the  rock  bridge  of  the  edge 
models  were  similar  to  those  obtained  from  the  full  model.)  To 
obtain  the  parameters  Eq  and  A,  we  went  through  several  calibration 
cycles  using  the  30°-75°  fracture  geometry.  In  other  words,  we 
estimated  Eo  and  A,  calculated  the  corresponding  parameter  Dsat  and 
used  the  finite  element  model  with  these  parameters  to  compute  the 
applied  load  at  which  DSat  is  reached  within  the  entire  rock  bridge. 
This  computed  load  was  compared  to  the  measured  load.  This 
calibration  led  us  to  chose  e0  =  0.0024  and  A  =  0.785  corresponding 
to  a  DSat  =  0.19.  We  then  used  the  same  parameter  values  for  other 
geometries.  The  computed  and  experimentally  obtained  loads  for  all 
non-overlapping  geometries  are  shown  in  Fig.  2.15.  As  this  figure 
shows,  there  is  a  satisfactory  agreement  between  the  experimental 
and  computed  values  with  discrepancies  that  are  roughly  the  same 
for  all  geometries. 

2.4  Concluding  Comments  and  Outlook 

The  research  on  fractures  of  fractured  rock  revealed  that  pre¬ 
existing  fractures  coalesce  in  two  different  modes  depending  on  their 
geometry.  If  the  pre-v  voting  fractures  overlap,  they  coalescence 
through  interconnection  of  the  developing  wing  cracks;  if  the  pre¬ 
existing  fractures  do  not  overlap,  coalescence  occurs  through 
secondary  cracks  which  occur  in  addition  to  and  after  the  wing 
cracks.  Such  fracture  coalescence  has  not  been  observed  so  far,  and, 
consequently,  no  analytical  prediction  model  exists.  (As  a  matter  of 
fact,  fracture  coalescence  analysis  for  overlapping  pre-existing 
fractures  only  predicts  wing  crack  geometry  but  not  the  actual 
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coalescence.)  The  coalescence  of  non-overlapping  fractures  occurs  at 
loads  which  are  lower  than  those  for  overlapping  fractures  with 
otherwise  identical  geometries.  For  this  reason,  we  developed  a  new 
analytical  model  which  we  validated  successfully. 

These  experimental  and  analytical  results  represent  thus  a 
significant  breakthrough  which  is  not  only  important  for  fractured 
rock  but  for  any  brittle  material  with  pre-existing  discontinuities 
such  as  cracks  or  fractures.  The  results  are  potentially  very  relevant 
to  the  representation  of  behavior  varying  from  microcrack 
coalescence  in  brittle  materials  to  fracture  and  fault  coalescence  in 
geology. 

Clearly,  only  a  first  step  has  been  made  in  this  research. 
Numerous  issues  remain  to  be  resolved.  Experimentally  it  is 
necessary  to  consider  other  materials  and  other  fracture  geometries. 
Eventually  we  will  also  have  to  consider  more  complex  loading 
conditions  and  the  possibility  that  different  types  of  fracture 
coalescence  occur  simultaneously  or  in  interrelated  sequence. 
Experiments  would  probably  gain  from  whole  field  measurement 
technology  such  as  holography  and  speckle  interferometry. 

Analytically  one  would  like  to  eventually  model  the  process  in 
more  detail.  In  particular,  it  would  be  desirable  to  represent  crack 
propagation  considering  the  combination  of  tensile  fractures  and 
shear  fracturing  that  evidently  takes  place. 


29 


3.  Centrifuge  Modelling 

3.1  Problem  Statement  and  Research  Schedule 

The  great  advantage  of  centrifuge  modelling  is  the  possibility 
to  apply  stresses  at  levels  corresponding  to  reality,  and  which  are 
geostatically  distributed  in  the  model.  In  preparation  for  the  present 
research,  we  have  conducted  an  extensive  study  of  scaled  modelling 
and  of  centrifuge  testing  in  rock  mechanics  [Literature  Review  of 
Geotechnical  Centrifuge  Modelling  with  Particular  Emphasis  on  Rock 
Mechanics,  1987  by  Joseph,  Einstein,  Whitman  and  Rock  Modelling 
Using  the  Centrifuge,  1987  by  Joseph  and  Einstein].  These  studies 
have  confirmed  the  questions  raised  about  the  applicability  of 
centrifuge  testing  in  modelling  of  jointed  rock,  i.e.  in  modelling  of 
discontinua.  The  questions  arise  from  the  possibility  that 
discontinuum  behavior  may  depend  on  absolute  displacements. 

Since  the  displacements  in  the  centrifuge  are  scaled  at  the  ratio  of 
model  to  prototype  geometry,  they  are  smaller  than  the  prototype 
displacements.  If  absolute  displacements  govern  the  behavior,  a 
centrifuge  test  will  therefore,  not  correctly  model  the  real  behavior, 
if  the  same  material  is  used  for  model  and  prototype. 

One  of  our  objectives  in  this  research  was  thus  to  determine 
the  applicability  of  centrifuge  modelling  for  discontinua  through  a  so 
called  "modelling  the  model"  approach.  In  this  approach  the 
behavior  of  differently  scaled  models  of  the  same  prototype  is 
compared.  If  absolute  displacements  (or  other  absolute 
characteristics)  rather  than  scaled  displacements  govern  the 
behavior,  this  will  become  apparent.  The  problem  we  chose  to  study 
is  arching  behavior  of  a  jointed  rock  mass  using  the  classic  trapdoor 
experiment. 

An  obvious  starting  point  for  trapdoor  arching  tests  with 
jointed  material  were  analogous  tests  with  granular  material.  This 
would  allow  us  to  compare,  and  thus  in  a  way  calibrate  our  trapdoor 
experimentation  with  the  many  trapdoor  experiments  on  granular 
material  which  have  been  performed  under  normal  gravity.  It  also 
provided  a  basis  for  comparison  with  the  few  centrifuge  trapdoor 
experiments  on  granular  material  conducted  by  others.  Differences 
between  trapdoor  1-g  and  elevated  g  experiments  with  such  material 
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would  have  significant  implications  on  analysis  and  design 
approaches  which  are  based  on  arching.  Following  the  trapdoor 
studies  on  sand  were  those  involving  a  discontiuuum  i.e  a  jointed 
rock  model.  As  mentioned  before  this  would  allow  us  to  investigate 
the  applicability  of  centrifuge  testing  to  model  discontinua.  It  would 
also  provide  us  with  better  knowledge  on  arching  behavior  in  such 
materials.  As  a  final  problem  area  in  this  research,  we  studied  the 
arching  concept  and  related  theories.  The  intention  was  to  reconcile 
our  experimental  observations  with  the  existing  theories  and,  if 
necessary,  develop  a  new  or  modified  theory. 

The  first  year  of  this  research  project  was  devoted  to  the 
development  of  a  trapdoor  setup  in  the  centrifuge  and  to  conducting 
initial  tests  with  sand  (see  "Stochastic  and  Centrifuge  Modelling  of 
Jointed  Rock,  Annual  Report,  1987/88).  In  the  second  year  of 
research  (see  "Stochastic  and  Centrifuge  Modelling  of  Jointed  Rock", 
Annual  Report  1988/89)  we  modified  the  experimental  setup  based 
on  our  initial  sand  experiments  and  conducted  the  remainder  of  the 
sand  experiments.  We  also  started  the  tests  on  discontinua  (jointed 
rock)  which  was  also  the  main  topic  of  the  third  year  research. 

During  this  third  year  we  also  conducted  tests  on  a  scaled  granular 
substance,  namely,  glass  beads  of  different  sizes.  During  the  second 
and  third  year,  parallel  research  at  the  Ruhr  University,  Bochum,  was 
performed  on  jointed  rock  with  a  somewhat  larger  trapdoor  and  joint 
spacing.  (This  research  which  was  funded  by  German  sponsors, 
enhanced  our  modelling  the  model  approach  and  led  to  a  paper, 
"Trapdoor  Experiments  With  Simulated  Jointed  Rock",  Iglesia,  et  al., 
by  both  research  groups.) 

3.2  Experimental  Arrangement 
3.2.1  MIT  Setup  and  Calibration 

The  essential  features  of  the  experimental  apparatus  at  MIT  are 
presented  in  Figure  3.1.  The  movement  of  the  trapdoor  was  achieved 
by  means  of  a  system  of  wedges.  The  lower  wedge  was  pulled  by  a 
motor-driven  rod  through  a  gear  box  assembly,  causing  the  upper 
wedge,  which  is  restrained  laterally,  to  slide  down  at  a  vertical 
displacement  rate  of  0.018  mm/sec.  The  vertical  force  on  the 


Figure  3.1:  MIT  Trapdoor  Apparatus 
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trapdoor  was  measured  by  load  cells  installed  between  the  trapdoor 
and  the  upper  weoge.  The  displacement  of  the  door  was  monitored 
by  a  transducer,  the  probe  of  which  was  attached  to  the  lower  wedge 
through  a  tiny  piece  which  protrudes  vertically  at  the  toe  of  the 
wedge.  (Note  that  the  measured  displacement  was  the  horizontal 
movement  of  the  lower  wedge;  the  vertical  translation  of  the 
trapdoor  can  be  computed  from  simple  geometric  relations.) 

The  trapdoor  strip  consisted  of  three  aluminum  segments,  each 
6.35  mm  thick  and  supported  by  a  ball  and  a  roller,  as  shown  in 
Figure  3.2.  One  segment,  4  inches  (101.6  mm)  in  length,  has  been 
placed  in  between  two  shorter  segments,  each  3  inches  (76.2  mm) 
long.  The  balls  transmitted  part  of  the  door  load  onto  the  load  cells 
through  cylindrical  steel  caps  which  were  positioned  right  at  the 
center  of  each  load  cell.  (To  ensure  a  truly  concentric  loading, 
another  set  of  balls  and  steel  caps  were  installed  underneath  the 
center  of  each  load  cell.)  The  roller  support  was  located  almost  at  the 
end  of  each  segment,  while  the  ball  support  (atop  a  load  cell)  was 
situated  at  a  third  of  the  segment  span  from  the  other  end.  The 
edges  of  the  trapdoor  segments  have  been  tapered  slightly  so  as  to 
minimize  friction  with  the  side  walls  of  the  trench.  In  the  absence  of 
other  measurements  related  to  the  door  load,  a  critical  assumption 
being  made  in  the  reduction  of  the  data  was  that  the  resultant  force 
due  to  the  earth  pressure  acts  at  the  center  of  each  door  segment. 

A  number  of  calibration  experiments  were  run  with  this  setup 
consisting  of  runs  with  no  load  and  with  calibrated  masses  on  the 
trapdoor.  After  initial  tests  with  sand  indicated  that  the  stiffness  of 
the  door  may  be  different  from  that  of  the  adjacent  base,  a  special 
series  of  tests  was  run  to  investigate  this.  As  it  turned  out,  the  door 
was  stiffer  than  the  adjacent  base  and  "attracts"  thus  more  load 
before  the  trapdoor  is  lowered.  This  explains  why  the  load 
displacement  curves  in  our  tests  started  at  loads  which  are  greater 
than  the  theoretical  load  Ngph  (Ng  =  acceleration  in  multiples  of 
gravity,  p=density,  h=overburden  depth).  A  correction  was  easily 
possible  by  shifting  the  load  displacement  curve  such  that  the 
theoretical  load  occurs  at  zero  displacement,  as  illustrated  in  Figure 
3.3.  The  shifting  of  the  curve  is  apparent  in  Figure  3.3a  and  3.3b. 
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The  portion  of  the  curve  where  the  ratio  (P/Po)  of  the  measured  load 
to  the  theoretical  load  is  less  than  or  equal  to  unity  (Figure  3.3c)  was 
then  considered  to  be  the  correct  curve  in  the  active  arching  mode. 

3.2.2  Ruhr  University  Bochum  (RUB')  Setup 

The  trapdoor  system  at  RUB,  depicted  in  Figure  3.4  was  slightly 
different  from  the  MIT  arrangement.  Movement  was  controlled  by  a 
hydraulic  system  connected  to  a  piston  which  supported  the 
trapdoor.  To  prevent  premature  piston  settlement  at  the 
acceleration  was  increased  prior  to  door  lowering,  back  pressure  was 
applied  to  the  piston  through  a  column  of  water  contained  in  a 
medium-high  pressure  vessel.  Before  dropping  the  piston,  the  back 
pressure  was  cut  off  by  closing  the  solenoid  valve  1,  and  the  piston 
was  subsequently  lowered  by  opening  valve  2. 

The  front  of  the  strong  box  was  formed  by  a  70  mm  thick 
perspex  face  through  which  deformations  of  the  model  were 
recorded  on  video  cassette.  In  the  RUB  test  package,  there  were  four 
trapdoor  segments,  each  10  mm  thick,  97.5  mm  long,  and  also  made 
of  aluminum.  The  segments  were  designated  number  1  through  4 
from  the  front  face  to  the  back  wall  of  the  strong  box.  Only  three  of 
these  -  segments  1,  2,  and  4  -  were  equipped  with  load  cells,  because 
segment  3  was  instrumented  with  a  displacement  transducer  to 
monitor  the  trapdoor  movement.  Each  segment  was  simply 
supported  as  in  the  MIT  setup,  except  that  a  "knife  edge",  was  used 
instead  of  a  roller. 

3.3  Tests  on  Granular  Material 

3.3.1  Tests  with  Sand 

For  these  experiments,  New  Jersey  4/14*  (coarse)  sand  with 
subangular  grains  (specific  gravity  Gs=2.66)  having  an  average  size  of 

2.1  mm  was  used.  The  typical  testing  procedure  as  shown  in  Figure 
3.5  consisted  of  applying  the  elevated  gravitational  level  in  four 
stages.  Once  the  highest  acceleration  level  was  reached,  the  trapdoor 

*  These  numbers  indicate  the  U.S.  Standard  sieve  sizes  through  which  the  sand 
particles  pass  and  on  which  they  are  retained.  Thus,  the  sand  used  passes 
through  the  #4  (4.76mm)  sieve  and  is  retained  on  the  #14  (1.41mm)  sieve. 
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was  moved  down  up  to  the  maximum  possible  displacement  of  4 
mm.  After  this,  the  acceleration  was  reduced  until  the  lg  level  was 
reached. 

The  results  showed  the  expected  decrease  of  the  trapdoor  load 
as  displacement  increased.  (Fig.  3.6  shows  typical  load  displacement 
curve  and  it  also  indicates  that  the  minimum  load  increased  as  the 
overburden  depth  decreased.  This  is  entirely  analagous  to  trapdoor 
experiments  under  normal  gravity  run  by  others  (e.g.  Terzaghi, 

1936).  The  normalized  displacement  (i.e.  displacement/trapdoor 
width)  at  which  the  minimum  load  occurred  was  about  0.02  for  the 
centrifuge  experiments  and  previous  tests  run  under  normal  gravity 
by  others.  What  was  different  in  high  g  centrifuge  testing  was  the 
fact  that,  after  reaching  the  minimum  load,  the  load  increase  with 
displacement  was  much  less  than  in  earlier  normal  gravity  tests. 

Modelling  the  model  was  investigated  by  properly  scaling 
geometry  and  gravity  level.  As  shown  in  Fig.  3.7  the  corresponding 
curves,  e.g.  2"  overburden  on  2"  door  at  40g  compared  to  1" 
overburden  on  1"  door  at  80  g  coincide. 

3.3.2  Tests  with  Glass  Beads 

Four  different  size  glass  beads  with  mean  diameters  of  0.5  mm, 
1.5  mm,  3.0  mm  and  6.0  mm  were  used.  They  were  placed  in  the 
same  manner  as  the  sand,  and  the  spinup,  door  lowering  and 
termination  of  the  centrifuge  tests  was  identical  to  those  on  sard. 
When  everything  was  scaled  appropriately,  the  normalized 
force/normalized  displacement  curves  in  Fig.  3.8  coincided.  This 
indicated  that  proper  scaling  should  include  the  scaling  of  grain  size. 
Subsequent  tests  were  run  in  which  either  grain  size  was  varied, 
while  the  trapdoor  width  was  held  constant  and  vice  versa  (door 
width  varied,  grain  size  constant);  these  tests  indicated,  however, 
that  the  scaling  of  the  door  size,  i.e.  of  the  "structure"  is  much  more 
important  than  the  scaling  of  the  grain  size. 
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3.3.3  General  Observations  on  Centrifuge  Tests  with  Granular 

Materials 

As  can  be  seen  in  the  two  preceding  sections,  the  load 
displacement  curves  first  dropped  to  a  minimum  load  followed  by  a 
slowly  increasing  load,  for  both  sand  and  glass  beads.  As  a  matter  of 
fact,  the  displacement  when  the  minimum  load  is  reached  was  1  to 
2%  of  the  door  width  and  was  approximately  the  same  in  all  tests  run 
by  us  in  the  centrifuge  at  various  g  levels  and  by  others  at  normal 
gravity.  What  was  different  in  the  high  gravity  centrifuge 
experiments  was  the  slower  load  increase  after  the  minimum  load 
compared  to  lg  experiments.  Also  significant  was  the  fact  that  an 
actual  arch  seemed  to  form,  this  arch  was  initially  curved  ("gothic" 
arch)  but  became  triangular  as  the  displacements  increased. 

3.4  Tests  with  a  Jointed  Medium 
3.4.1  Test  Layout 

Initial  tests  were  run  with  balsa  wood  rods  of  square  cross- 
section  and  dimensions  of  1/4"  and  1/2"  square.  The  results 
scattered  significantly  and  made  these  tests  difficult  to  interpret. 
Nevertheless,  arching  effects  occurred  leading  to  a  load  reduction 
down  to  25%  of  the  original  load. 

The  main  test  series  on  discontinua  both  at  MIT  and  at  RUB 
were  run  with  aluminum  rods.  Ideally,  the  rods  should  be  made  of  a 
rock-like  material,  but  it  is  difficult  to  find  one  which  can  easily  be 
cast  or  cut  into  small  pieces  of  uniform  size.  Aluminum  rods  were 
used  since  they  are  readily  available  in  sizes  small  enough  such  that 
a  reasonable  number  of  pieces  fit  on  the  relatively  narrow  trapdoor 
strip.  Also,  if  the  main  concern  is  the  particular  displacement 
dependent  shearing  between  planar  surfaces  of  semi-rigid  elements 
and  not  the  fracturing  phenomena,  experiments  with  metal  rods 
should  provide  meaningful  results.  Moreover,  aluminum  is  well 
suited  for  simulating  rock  due  to  its  high  strength-to-weight  ratio 
(ideal  for  modelling  semi-rigid  components)  and  a  specific  gravity  of 
2.66  which  is  not  too  different  from  that  of  rock. 

The  jointed  rock  models  were  composed  of  rods,  square  in 
cross-section,  which  were  juxtaposed  side  by  side  and  top  of  each 
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other  to  simulate  joint  sets  withuniform  spacing.  Two  types  of  stack- 
up  schemes  were  implemented:  one  in  which  the  rods  were  directly 
placed  on  top  of  one  another,  and  the  other  in  a  staggered  fashion 
like  bricks  (Figure  3.9).  Following  a  "modelling-of-models" 
technique,  different  but  geometrically  similar  sizes  of  small-scale 
models  (i.e.,  aluminum  rods)  and  of  trapdoors  were  used  with 
corresponding  accelerations  in  the  centrifuge.  At  MIT,  the  door 
widths  were  2"  (50.8  mm),  1"  (25.4mm)  and  1/2"  (12.7mm)  for  the 
1/2"  (12.7mm)  [  ]  l/4"(6.35mm)[  ],  and  1/8"  (3.18  mm)[  ]  rods, 
respectively,  and  the  tests  were  carried  out  at  40g,  80g,  and  160g. 

The  RUB  tests  are  performed  at  25g  and  50g  using  corresponding 
door  widths  of  80mm  and  40  mm  for  aluminum  rods  sizes  of 
20  mm  [  ]  and  10mm[  ],  respectively.  (Since  the  experimental  setup 
in  this  study  was  presumed  to  be  in  plane  strain,  the  dimensions 
along  the  length  of  the  trapdoor  were  considered  to  be  irrelevant  and 
were,  thus,  not  scaled.) 

The  lateral  confinement  of  the  aluminum  rods  was  also 
standardized  for  both  the  MIT  and  RUB  tests.  Specifically,  silica 
quartz  glass  beads,  fairly  uniform  in  size  and  with  a  specific  gravity 
of  2.55,  were  placed  on  the  sides.  The  use  of  glass  spheres  also 
allowed  one  to  vary  the  diameter  of  the  beads  in  accordance  with  the 
scaling  of  the  rod  and  door  sizes.  The  three  glass  bead  diameters 
utilized  at  MIT  were  6  mm,  3  mm,  and  1.5  mm  for  the  40g,  80g,  and 
160g  tests,  respectively.  At  RUB,  5  mm  diameter  glass  beads  were 
used  for  the  50g  experiments  and,  due  to  the  unavailability  of  10 
mm  diameter  beads,  for  the  25g  tests  as  well.  A  listing  of  the 
pertinent  series  of  tests  conducted  at  MIT  and  at  RUB  is  given  in 
Table  3.1. 

3.4.2  Direct  Stack  Tests 

The  results  from  the  trapdoor  experiments  with  the  direct 
stack  layout  are  displayed  in  Figures  3.10  and  3.11,  The  final 
configuration  -  the  appearance  of  the  blocky  system  at  the  end  of 
each  test  resembles  Figure  3.10.  That  is,  the  four  columns  of  rods 
directly  above  the  trapdoor  simply  followed  the  vertical  door 
movement.  In  Figure  3.11,  the  measured  force  on  the  trapdoor, 


Rod  Size 

Total  Height 

Door  Width 

Glass  Bead  Size 

Gravity  Level 

MIT 

1. 

\  in 

5.00  in 

2.0  in 

6.0  mm 

40g 

(12.7  mm) 

(127  mm) 

(50.8  mm) 

2. 

}in 

2.50  in 

1.0  in 

3.0  mm 

80g 

(6.35  mm) 

(63.5  mm) 

(25.4  mm) 

3. 

S  in 

1.25  in 

0.5  in 

1.5  mm 

160g 

(3.18  mm) 

(31.8  mm) 

(x2.7  mm) 

RUB 

4. 

10  mm 

100  mm 

40  mm 

5.0  mm 

50g 

5. 

20  mm 

200  mm 

80  mm 

5.0  mm 

25g 

Table  3.1.  Testing  Program 
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b.  Brick  Stack 


Figure  3.9:  Jointed  Medium  Stack-up  Patterns 


Figure  3.10:  Direct  Stack  Tests:  Final  Configuration 
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M  I  T  experiments: 
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Figure  3.11:  Direct  Stack  Tests:  Normalized  Load  vs.  Displacement  Data 
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normalized  with  respect  to  the  weight  (at  Ng)  of  the  material  directly 
above  the  door,  has  been  plotted  against  the  vertical  door  translation 
for  the  various  tests. 

Note,  in  the  MIT  experiments,  that  the  measured  force  on  the 
door  at  zero  displacement  was  greater  than  the  theoretical  geostatic 
force.  This  was  caused  by  the  greater  stiffness  of  the  trapdoor 
relative  to  the  adjacent  base  as  discussed  in  Section  3.2.1.  The 
results  show,  with  minor  experimental  scatter,  that  upon  lowering  of 
the  door,  an  immediate  reduction  in  load  ensued.  This  was  followed 
by  a  stage  of  approximately  constant  load  up  to  a  door  displacement 
of  1  to  1.25mm,  then  a  further  decrease  in  load  to.  about  1.5  mm 
displacement,  and  a  gradual  force  increase  toward  another  fairly 
constant  value. 

3.4.3  Brick  Stack  Test 

For  the  brick  stack  cases,  the  results  are  given  in  Figures  3.12 
and  3.13.  In  the  final  configuration  (Figure  3.12),  a  triangular 
pattern  is  typically  formed  in  the  bottom  four  layers  just  above  the 
trapdoor.  This  is  reflected  in  the  measured  loads  (Figure  3.13)  which 
became  constant  towards  the  end  of  each  experiment.  The 
normalization  of  the  forces  in  Figure  3.13  was  done  with  respect  to 
the  theoretical  geostatic  weight  of  a  volume  of  material  above  the 
door  bounded  by  vertical  planes  passing  through  the  edges  of  the 
door.  A  maximum  of  40  blocks  could  be  placed  above  the  door,  and 
the  10  rods  that  followed  the  door  movement  in  each  test  constitute 
thus  25  percent  of  the  theoretical  geostatic  force. 

As  the  trapdoor  settled  in  the  brick  stack  case,  the  transition 
from  the  initial  force  to  the  minimum  load  followed  a  much  smoother 
curve  than  in  the  direct  stack  experiment.  It  is  evident  from  Figure 
3.13  that  the  minimum  load  occurred  at  about  the  same  amount  of 
door  displacement  for  all  geometrically  similar  tests.  This 
observation  is  even  more  pronounced  if  the  data  are  adjusted  such 
that  zero  displacement  is  taken  where  the  normalized  force  ratio 
equals  unity,  shifting  the  curves  in  Figure  3.13  horizontally  to  those 
in  Figure  3.14. 
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Figure  3.12:  Brick  Stack  Tests:  Final  Configuration 
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M  I  T  experiments:  * 

rdls,  6”  high,  on  2’’  door  •?  iOg 
0  -  J"c  rods,  2.5"  high,  on  1"  door  §  80g 
X  -  rods,  1.25”  high,  on  i”  door  160g 

RUB  experiments: 

□  -  10  lmnD  rods,  100  inm  high,  on  40  nun  door  ■§» 

0-20  minQ  rods,  200  mm  high,  on  80  mm  door  ?' 
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50g  4 
25g  4 
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Figure  3.14:  Brick  Stack  Tests:  Normalized  Load  vs.  Adjusted  Displacement 
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3.4.4  Interpretation  of  Tests  Results 

The  results  from  the  "modelling-of-models"  exercise  with 
discontinua  showed  that  the  geometrically  similar  setups  at 
appropriately  scaled  g-levels  exhibited  similar  load-displacement 
behavior  in  the  trapdoor  tests.  The  force  values,  especially  the 
minimum  loads,  conformed  approximately  with  the  usual  centrifuge 
scaling  laws.  However,  these  forces  occurred  at  the  same  amount  of 
absolute  door  displacement,  in  violation  of  the  same  scaling  rules. 
Similitude  requirements  for  scaling  of  length  dimensions  are  simply 
not  met  by  the  observed  displacements. 

The  fact  that  shearing  depends  on  the  absolute  displacement 
between  surfaces  provides  a  serious  obstacle  to  centrifuge 
experimentation  on  models  involving  discontinuities.  It  may  be 
extremely  difficult,  if  at  all  possible,  to  extrapolate  results  from 
model  tests  to  actual  prototype  conditions  unless  the  displacement 
dependence  is  also  considered,  possibly  by  scaling  roughness  on  the 
discontinuity  surfaces. 

3.5  Arching  Theories 

The  settling  of  a  trapdoor  underneath  a  geomaterial  produced 
an  arch  just  above  the  door.  The  slope  of  this  arch  is  initially 
parabolic  Fig.  3.15a,  then  becomes  triangular  (3.15b).  Eventually,  at 
large  trapdoor  displacement,  the  arch  collapses  and  a  prismatic  body 
(Fig.  3.16)  of  the  overlying  material  folllows  the  movement  of  the 
trapdoor.  In  our  centrifuge  experiments,  which  are  limited  to  small 
displacements,  we  were  able  to  observe  the  initially  curved 
(parabolic)  arch  in  both  the  sand  and  glass  bead  tests,  as  well  as  a 
transition  to  a  triangular  arch  in  the  glass  bead  tests.  Prismatic 
bodies  at  very  large  trapdoor  displacements  in  centrifuge  tests  were 
observed  by  Stone  (1988). 

This  mechanism  of  arching  is  consistent  with  well  established 
relationships  between  (incremental)  shearing  strains  and 
(incremental)  volumetric  strains  in  geomaterials.  At  small  strains,  a 
granular  medium  tends  to  dilate  or  increase  in  volume  depending  on 
the  stress  level  and  porosity.  This  ability  of  granular  soils  to  carry 
shearing  stresses  brings  about  the  formation  of  the  arch  above  the 


(b)  Triangular  Arch 


Figure  3.15:  Arching  at  Relatively  Small  Displacements 
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trapdoor,  and  considerable  dilation  of  the  soil  occurs  soon  after  the 
door  is  Hi  solaced  At  large  strains,  the  geomaterial  tends  to  deform 
with  haiw  any  change  in  volume.  The  eventual  collapse  of  the  arch 
leading  to  vertical  failure  surfaces  is  a  manifestation  of  this  stage  of 
isochoric,  or  constant-volume,  deformation. 

The  formation  of  an  arch  above  a  trapdoor  or  more 
importantly,  above  the  crown  of  a  tunnel  has  been  observed  or  at 
least  assumed  in  a  number  of  analysis  methods  for  tunnel  support. 
Engesser  (1881)  and  Bierbaumer  (1913)  assumed  parabolic  or 
triangular  arches,  respectively  (Fig.  3.15)  while  Terzaghi  worked 
with  the  ultimate  prismatic  body  (Fig.  3.16).  Evans  (1983)  developed 
an  arching  analysis  based  on  observed  triangular  arching  and 
rigorous  consideration  of  plasticity  theory  (Volume  II  of  this  report 
series  contains  a  detailed  review  of  these  and  other  theories). 
Predictions  with  these  methods  are  compared  in  Table  3.2  (taking 
<j)=39°  as  determined  for  the  sand  used  on  our  centrifuge  tests). 


H/B 

Engesser  K=Ka 

Terzaghi  K=1 

Evans 

(P/P)min2 

(K=l  2) 

(P/Po-lult 

1 

P/Po 

P/Po 

0.309 

0.516 

0.258 

0.495 

2 

0.133 

0.297 

0.154 

0.315 

OH 

0.068 

0.168 

0.077 

0.182 

Table  3.2.  Comparisons  of  different  theories  applied  to  the 

centrifuge  tests  with  sand:  Computed  values  of  force  on 
trapdoor  normalized  with  respect  to  overburden  (P/P0). 
(K=Lateral  earth  pressure  coeff.;  Ka=  active  earth  pressure 
coeff.;  H/B=ratio  of  overburden  height  to  trapdoor  width) 


In  our  tests  the  mimimum  normalized  P/P0  forces  were  0.36,  0.18, 
0.09  for  H/B=l,2,4  respectively.  Terzaghi's  theory  which  is  based  on 
the  ultimate  load  naturally  produces  larger  loads.  Of  the  other  two, 
Evans  predictions  seem  to  be  closer  to  our  observations. 
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Consequently,  we  formulated  a  theory  by  combining  the 
formulations  of  Engesser  (1881),  Bierbaumer  (1913)  and  Evans 
(1983).  The  proposed  expression  for  the  minimum  load  is 


2  K*  cot  (b 

>  .  =  yB  ( -  ■  »  -  +  — T1) 

"""  2cot*  +  |K  6 


(3.1) 


where  B  is  the  width  of  the  trapdoor,  H  is  the  depth  of  overburden,  y 
is  the  unit  weight  of  the  soil,  (j)  is  the  friction  angle,  and 

cos 2  6 

K>- - rr~ 

1  +  sin  z(j) 


The  beauty  of  this  new  theory  lies  in  the  fact  that  one  need  not 
assume  a  value  for  the  lateral  coefficient  of  pressure  K  -  the  main 
source  of  controversy  in  the  other  methods.  It  comes  about 
automatically  from  a  Mohr  circle  construction  by  hypothesizing  that 
the  arch  itself  is  a  failure  surface,  with  a  slope  from  the  vertical 
through  the  edges  of  the  trapdoor  equal  to  the  friction  angle  of  the 
soil.  The  shape  of  the  arch  is  assumed  to  be  a  parabola,  and  the  force 
on  the  trapdoor  is  the  sum  of  the  weight  of  material  below  this  arch 
(which  depends  on  B,  y  and  4>)  and  the  effect  of  the  induced  stresses 
from  the  sides  due  to  the  load  redistribution  (which  is  a  function  of 
H,  y  and  <(>). 

Computations  in  this  study  have  shown  that  the  predicted 
minimum  load  on  the  structure  (trapdoor)  is  approximately  the  same 
whether  the  shape  of  the  arch  is  curved,  as  in  Engesser's  theory,  or 
triangular,  as  in  Bierbaumer's  or  Evan's  theory.  Both  shapes  have 
indeed  been  observed  experimentally,  with  the  triangle-like  pattern 
occuring  at  larger  displacements.  Thus,  even  if  the  arch 
configuration  transforms  from  a  curved  to  a  triangular  shape  while 
the  door  is  moved  down,  the  load  on  the  door  is  apparently 
unchanged.  This  may  explain  why  the  measured  forces,  as  soon  as 
the  minimum  is  reached,  stay  at  more  or  less  the  same  level  for  a 
considerable  amount  of  displacement. 


The  trapdoor  experiments  both  with  granular  and  jointed 
media  have  shown  some  interesting  similarities.  Specifically,  the 
formation  of  a  physical  arch  in  the  former  is  analogous  to  the  brick 
stack  case.  As  observed  by  Stone  (1988)  in  his  centrifuge  tests,  the 
rupture  surfaces  tend  to  become  vertical  as  the  trapdoor  is 
continually  lowered  -  just  like  in  the  direct  stack  arrangement. 

3.6  Concluding  Comments  and  Outline 

The  three  objetives  of  this  research,  to  study  arching  in 
granular  and  jointed  media  in  the  realistic  geostatic  stress 
environment  of  the  centrifuge,  to  examine  the  applicibility  of 
centrifuge  testing  through  a  modelling  the  model  approach  and  the 
development  -  if  necessary  -  of  a  new  arching  theory  were  achieved. 
Arching  leads  to  a  load  reduction  as  the  trapdoor  displaces.  In 
previous  arching  experiments  under  normal  gravity  a  minimum  load 
was  reached  followed  by  a  relatively  rapid  load  increase  afterwards 
while  in  the  centrifuge  at  high  g  levels  the  load  increased  very 
slowly  after  the  minimum.  This  can  have  a  significant  effect  on 
tunnel  design  where  one  assumes  a  load  reduction  with  displacement 
but,  so  far,  has  been  very  conservative  as  far  as  the  admissable 
displacements  are  concerned.  The  reason  for  this  conservativism  is 
the  possibility  of  loosening,  i.e.  the  development  of  the  higher 
ultimate  load.  As  shown,  this  ultimate  load  develops  very  slowly 
under  realistic  stress  conditions. 

We  were  able  to  pretty  conclusively  determine  that  centrifuge 
testing  is  applicable  in  modelling  with  granular  materials  while  the 
applicability  is  limited  in  jointed  media.  The  fact  that  shearing 
resistance  in  discontinua  depends  on  absolute  and  not  relative 
displacements  is  the  reason  for  this. 

The  trapdoor  tests  and  a  number  pre-existing  arching  theories 
allowed  us  to  develop  a  new  arching  theory  which  considers  the  fact 
that  a  real  arch  of  parabolic  or  triangular  shape  occurs  above  the 
trapdoor.  The  load  is  independent  of  the  lateral  stress  coefficient  (in 
contrast  to  other  arching  theories). 
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Future  work  will  have  to  extend  this  arching  theory  to  the 
ultimate  load.  In  this  context,  physical  experiments  in  the  centrifuge 
will  be  needed  to  accurately  establish  the  transformation  of  the  arch 
into  a  prismatic  body.  Inspite  of  the  fact  that  the  applicability  of 
centrifuge  testing  to  discontinua  is  limited,  cases  in  which  the 
geostatic  stress  conditions  may  be  more  important  than  displacement 
related  phenomena  need  to  be  studied.  Examples  are  rigid 
block/wedge  concepts  as  used  in  slope  stability  analysis. 
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4.  Stochastic  and  Topological  Fracture  Geometry  Model 
4.1  Problem  Statement  and  Possible  Approaches 

Attempts  at  describing  the  geometric  patterns  of  rock  joints  are 
numerous.  During  the  past  10  to  15  years,  much  of  this  effort  has 
been  directed  toward  the  development  of  stochastic  models,  and  the 
MIT  rock  mechanics  group  has  played  a  leading  role.  The  basic 
concept  of  these  models  is  to  statistically  sample  geometric  joint 
characteristics,  such  as  orientation,  trace  length  and  spacing 
(location)  followed  by  a  probabilistic  model  representation  of  the 
sampled  data.  Ideally  the  modelled  and  real  pattern  should  coincide 
(within  the  acceptable  statistical  range).  Given  the  uncertainty  of 
joint  geometry  (variability  or  range  of  characteristics;  limited  access 
through  exploration)  probabilistic  or  stochastic  approaches  represent 
in  most  cases  the  only  possible  way  to  come  up  with  engineering 
solutions  to  flow,  -stability,  -  and  deformability  problems.  Many 
geological  problems  also  rely  on  stochastic  representation  of 
geometry  for  their  solution. 

However,  the  existing  models  have  to  resort  to  significant 
simplifications,  specifically: 

•  They  (with  some  exceptions)  do  not  account  for  spatial 
nonhomogeneities  such  as  fracture  clustering. 

•  The  models  are  only  loosely  tied  to  the  geologic  genesis  of 
the  fractures.  In  particular,  most  models  assume 
independence  among  fracture  sets.  From  a  physical 
viewpoint,  this  assumption  is  often  incorrect. 

•  Only  in  a  few  cases  have  the  models  been  validated  using 
actual  fracture  data. 

As  a  consequence,  some  modelled  (simulated)  joint  patterns 
closely  represent  the  real  conditions  while  others  do  not.  It  is  also 
desirable  to  model  joint  geometry  in  a  way  that  corresponds  to  joint 
genesis.  Individual  joints  and  joint  sets  grow  in  a  mechanically 
based  sequence.  The  geometric  model  should  be  capable  of 
producing  joints  and  joint  sets  in  such  a  sequence.  Eventually,  one 
will  then  be  able  to  link  the  geometric  and  mechanical  models  to 
represent  joint  patterns  at  particular  sites. 
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What  we  were  trying  to  do  in  this  research  is  to  develop 
stochastic  models  which  more  realistically  represent  joint  patterns 
than  the  existing  ones,  and  which  allow  one  to  model  the  sequence  of 
joint  genesis.  (Eventually,  such  geometric  models  will  be  combined 
with  mechanical  models  such  as  those  obtained  from  the  research 
discussed  in  Section  2,  above). 

Figure  4.1  shows  typical  joint  trace  maps,  i.e.  the  intersections 
of  joints  with  the  outcrop  surface.  There  are  three  basic  approaches 
to  model  joint  patterns  such  as  the  ones  shown  in  Fig.  4.1. 

•  Point-based  processes.  For  instance  the  midpoint  distribution 
of  joint  traces  is  modelled.  It  is  then  possible  to  have  "marks" 
on  each  midpoint  which  describe  the  length  (size),  orientation 
and  other  characteristics  of  the  particular  joint. 

«  Fiber  processes.  The  trace  distribution  is  directly  modelled 
but  marks  can  again  be  added  to  characterize  other  properties 
such  as  size  in  the  third  dimension,  and  roughness. 

•  Random  closed  sets  (RACS)  model  the  joints  as  planar  features 
(i.e.  the  distribution  of  planar  features).  Again,  other 
properties  such  as  nonplanarity  can  be  associated  with  each 
feature. 

These  three  approaches  essentially  go  from  the  one 
dimensional  to  the  three  dimensional.  If  the  basic  approach  is  simple 
(one  dimensional)  other  characteristics  need  to  be  associated  through 
separate  processes  (The  "marks"  discussed  above).  So  far  we  have 
looked  into  and,  where  appropriate  developed,  point  and  fiber 
processes. 

In  our  investigation  of  point  processes,  we  tested  the  1) 
homogeneous,  isotropic  Poisson  point  process.  2)  inhomogeneous 
Poisson  point  process;  3)  Poisson  cluster  process;  4)  Cox  (doubly 
stochastic)  process;  5)  simple  inhibition  process;  6)  thinning  process; 
7)  soft  core  process. 

In  figure  4.2  the  midpoints  of  the  joint  traces  within  the 
rectangular  area  of  Fig.  4.1a  are  shown.  Tests  of  complete  spatial 
randomness  based  on  interevent  distances  and  nearest  neighbor 
distances  as  well  as  second  moment  properties  indicate  that  the 
points  in  Fig.  4.2  do  not  follow  a  homogeneous  Poisson  process.  The 


Figure  4.1a.  Map  of  Florence  Lake  Outcrop  (After  Segall  and 

Pollard,  1983;  Only  Traces  Within  Rectangular  Regi 
are  Considered) 
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Figure  4.1b.  Trace  Man  of  PA  100  ^rom  n-  -n,*  T  -no-, 

■-,w'  v*iv!in  uai  lOu  ana  parson,  1985) 


Figure  4.2 


Mid  Point  Representation  of  Traces  of  Florence  Lake 
Outcrop. 


64 


Inhomogeneous  Poisson  Point  Model  and  the  Doubly  Stochastic  (Cox) 
Model  were  applied.  With  the  inhomogeneous  Poisson  Point  Model 
we  also  included  the  joint  trace  length  as  a  mark  on  each  point. 
Nevertheless,  this  model  requires  that  joint  intensity  variations  be 
estimated  with  kernel  functions.  These  kernel  functions  are  not 
explicit  and  make  thus  estimation  of  joint  patterns  outside  the 
narrow  observed  area  difficult.  This  problem  docs  not  exist  with  the 
doubly  stochastic  process,  where  the  intensity  function  is  an  explicit 
stochastic  function.  Fig.  4.3  shows  the  midpoints  produced  with  the 
Cox  process  which  produces  a  pattern  very  similar  to  that  shown  in 
Fig.  4.2.  The  statistics,  such  as  the  second  moment  analysis,  confirm 
this. 

In  the  area  of  fiber  or  line  process  models  we  also  investigated 
a  number  of  possibilities  such  as  the  homogeneous  Poisson  fiber 
process,  the  parent  daughter  model,  the  fractal  model,  the  branching 
model,  the  hierarchical  model,  the  percolation  model  and  the  crack 
tessellation  model. 

4.2  Hierarchical  Fracture  Trace  Model 
4.2.1  Basic  Concept 

On  the  basis  of  our  research  into  point  and  fiber  processes,  we 
then  developed  a  hierarchical  fracture  trace  model  which  allows  one 
to  realistically  model  several  sets  of  fracture  traces  i.e.  the 
intersections  of  fractures  with  an  outcrop.  So  far,  we  have  developed 
the  model  to  handle  two  fracture  trace  sets,  but  an  expansion  to 
additional  sets  will  be  relatively  simple. 
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Figure  4.3. 
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Mid  Point  Pattern  Using  Cox  Process  for  Fracture 
Traces  Shown  in  Figure  4.2. 
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The  hierarchical  modelling  strategy  consists  of  the  following 

steps: 

1.  Model  the  midpoints  of  the  first  fracture  set  by  a  doubly 
stochastic  (Cox)  process,  which  allows  one  to  represent  spatial 
clustering. 

2.  Model  the  fracture  trace  length  and  orientation  of  the  first  set, 
including  spatial  correlation,  if  needed. 

3.  Test  for  dependence  between  the  location  (midpoints)  of  the 
second  fracture  set  and  the  traces  of  the  first  set.  If  there  is 
dependence,  account  for  it  when  modelling  the  location  of  the  second 
set;  otherwise,  model  the  second  set  independently,  according  to 
steps  1  and  2. 

4.  Model  trace  length  and  orientation  of  Set  2.  Determine  the 
probability  of  terminations  at  intersections  of  Set  2  fractures  with 
Set  1  fractures  and  model  the  terminations. 

5.  Repeat  for  subsequent  fracture  sets  (In  the  detailed 
development  and  example  application  included  in  this  paper,  only 
two  sets  will  be  considered). 

Steps  3,  4  and  5  reflect  the  hierarchical  aspect  of  the  model. 
The  example  used  for  illustration  is  pavement  PA  100  at  Yucca 
Mountain  (Fig.  4.1b)  .  A  straight-segment  idealization  of  the  trace 
pattern  is  shown  in  Fig.  4.4.  Two  characteristic  sets  can  be 
distinguished;  Set  1  includes  the  joints  with  nearly  parallel 
orientations  and  very  few  intersections.  Set  2  collects  the  remaining 
traces,  except  for  the  intensely  fractured  zones.  The  decomposition 
into  trace  sets  and  the  resulting  patterns  of  midpoints  are  shown  in 
Fig.  4.5. 
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Figure  4.5c.  Mid  Point  Map  of  PA100:  Mid  Points  of  Set  1  Traces  Only. 
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Figure  4.5d.  Mid  Point  Map  of  PA100:  Mid  Points  of  Set  2  Traces  Only 
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4.2.2.  Modelling  of  the  First  Trace  Set 

As  mentioned  before,  we  have  selected  a  Cox  point  process 
which  was  fitted  to  reproduce  second-moment  properties  of  the 
sample.  More  specifically  we  have  aimed  at  matching  the  sample 
value  of  the  function  Kjj(t).  The  indices  i  and  j  denote  in  general  two 
sets  of  points,  which  for  Kl  l  coincide  with  the  midpoints  of  Set  1 
traces.  For  an  unbounded  region  and  a  homogeneous  Poisson  process 
with  index  i,  Kii(t)  =  rct^,  and  for  any  two  independent  joint  processes 
i  and  j,  Kjj(t)  =  Tit2.  This  is  why  the  estimates  Kii(t)-7tt2  and  Kij(t)-7ct2 
can  be  used  to  detect  respectively  deviations  from  Poisson 
homogeneity  of  a  single  point  process  and  dependence  between  two 
point  patterns.  For  regions  with  complex  geometry,  the  expectations 
of  Kii(t)  -  jct^  and  Kij(t)-7tt^  under  the  above  conditions  are  not  zero 
and  are  best  obtained  through  Monte  Carlo  simulation.  The  Monte 
Carlo  approach  can  be  also  used  to  devise  acceptance  regions  at  any 
specified  significance  level.  For  example,  Figure  4.6  shows  the  mean 
value  and  the  acceptance  band  (at  significance  level  a=  0.02)  of 
Kn(t)-7ct2>  obtained  through  simulation  of  Set-1  midpoints  using  a 
specific  Cox  process.  Superposed  on  the  same  figure  is  the  empirical 
function  Kn(t)-7ct2,  obtained  from  the  data  in  Figure  4.5c.  The 
agreement  is  quite  good  over  the  entire  range  of  distances  at  which 
clustering  is  present.  For  comparison  with  Figure  4.5c,  Figure  4.7 
shows  a  realization  of  the  fitted  Cox  process  inside  the  actual  outcrop 
area. 


Next,  the  distribution  of  trace  length  is  considered.  As  is  well 
known  the  empirical  distribution  of  trace  length  for  the  scan  line 
method,  may  be  severely  biased.  Biases  exist  also  in  samples  from 
outcrops  such  as  that  of  Figure  4.1b.  A  general  method  to  correct  for 
these  biases  is  to  obtain  distribution  parameter  estimates  by  the 
method  of  maximum  likelihood. 

For  this  purpose,  one  must  classify  outcrop  traces  into  three 
types  T  and  measure  their  observed  length  L'  and  "window  of 
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Figure  4.6.  Second  Moment  Analysis  of  Set  1.  Acceptance  band 
(MIN  to  MAX)  and  average  value  (AVRG)  obtained 
with  model  simulations  are  shown.  The  second 
moment  mapped  data  (MAPPED)  is  also  given 
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observation"  w,  as  illustrated  in  Figure  4.8.  (Notice  that,  due  to 
censoring,  the  actual  length  L  of  type  2  and  type  3  traces  exceeds  the 
measured  length  L').  The  likelihood  of  the  trace  length  distribution 
Fl  can  be  expressed  as  a  function  of  the  sampled  values  Ti,Wi,LY,  and 
the  likelihood  can  be  maximized  to  obtain  estimates  of  the 
parameters.  For  example,  Figure  4.9  shows  the  logarithm  of  lengths 
7'  of  Set  1  as  a  function  of  the  mean  trace  length  m,  assuming  that 
the  Set  1  traces  have  exponential  length  distribution.  The  maximum 
likelihood  estimate  of  m  is  about  7.0  m.  It  is  interesting  to  notice 
that  consideration  of  only  Type  1  (fully  exposed)  traces  would  have 
given  a  much  smaller  estimate  (2.75  m);  see  Figure  4.9. 

For  Set  1,  we  have  not  modeled  the  variability  of  trace 
orientation,  since  this  variability  is  very  modest.  If  it  were 
important,  orientation  variability  would  be  modelled  analogously  to 
what  will  be  shown  below  for  Set  2. 

4.2.3  Modelling  of  the  Second  Trace  Set 

Considerable  complications  may  arise  in  the  modelling  of  trace 
Set  2,  due  to  dependence  on  the  fractures  of  Set  1.  One  should  start 
the  modeling  process  by  testing  independence  between  the  two  sets 
of  midpoints,  for  example,  using  the  function  K12  described  above. 
Whon  applied  to  the  data  in  Figs.  4.5c  and  4.5d,  the  test  indicates 
significant  interaction  and  more  specifically  a  tendency  of  Set  2 
points  to  be  located  near  Set  1  points.  From  a  physical  point  of  view, 
it  is,  however,  more  meaningful  to  consider  interaction  between  the 
location  (here,  the  midpoints)  of  Set  2  and  the  traces  (not  just  the 
midpoints)  of  Set  1.  For  this  purpose,  the  point  processes  of 
midpoints  of  Set  2  were  considered  with  intensity  functions  X,2(x,y) 
which  include  the  fact  that  the  midpoint  locations  of  Set  2  depend  on 
the  location  of  Set  1  traces.  K12  can  then  be  computed  using  such  a 
(theoretical)  intensity  function.  Figure  4.10a  shows  the  match 
between  the  theoretical  and  expirical  K12  functions.  The  fit  is 
obviously  good  and  the  empirical  function  remains  far  from  the 
bounds  of  the  acceptance  region  at  the  0.02  confidence  level. 


Figure  4.9.  Log-Likelihood  of  Mean  Trace  Length  for  Set  1, 
Assuming  Exponential  Distribution. 
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Figure  4.10a.  Bivariate  Second  Moment  Analysis  for  Sets  1  and  2  - 
Mid  Point  Correlation  based  on  the  Nearest  Neighbor 
Function  Method.  Acceptance  band  (MIN  to  MAX)  and 
Average  Value  (AVRG)  are  obtained  with  model 
simulations,  while  the  mapped  pattern  (MAPPEr  )  is 
from  actual  data. 
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igure  4.10b.  (Univariate)  Second  Moment  Analysis  for  Set  2  Mid 
Points  Using  the  Nearest  Neightbor  Function  Method 
Symbols  are  the  same  as  those  in  Figure  4.10a. 
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However,  the  model  still  needs  validation,  because  the  kJ2  function 

was  used  in  the  estimation  of  the  function  K.  A  specific  concern  is 
that  the  spatial  variation  of  X2  might  not  completely  explain  the 
clustered  pattern  of  Set  2  points.  A  verification  can  be  made  by 
considering  the  function  K22O):  this  function  is  quite  different  from 
K 12  and  therefore  is  an  appropriate  choice  for  validation.  A 
comparison  of  K22(t)-7tt2  with  results  from  simulating 
nonhomogeneous  Poisson  patterns  with  intensity  X-2  is  shown  in 
Figure  4.10b,  while  Figure  4.11  reproduces  a  typical  simulated 
pattern  which  should  be  compared  to  the  real  data  in  Fig.  4.5d. 

These  results  clearly  confirm  the  model. 

The  estimation  of  the  trace  length  distribution  of  the  second  set 
was  done  as  for  the  first  fracture  set.  The  result  is  that  such  a 
distribution  can  be  taken  as  exponential,  with  a  mean  value  around  2 
meters. 

Orientation  of  Set  2  traces  was  modelled  by  fitting  either  a  von 
Mises  or  a  wrapped  normal  distribution. 

Regarding  termination  of  Set  2  traces  on  Set  1  traces, 
T-intersections,  i.e.  terminations  and  X-ir.tersections,  i.e.  crossings  are 
considered.  The  empirical  probability  of  T  and  X  intersections 
(P[X]=1-P[T])  were  determined.  These  probabilities  have  then  been 
used  as  follows  (Fig.  4.12):  t2  is  a  generic  simulated  type-2  trace,  and 
A  and  B  are  intersection  points  of  t2  with  type-1  traces  closest  to  its 
endpoints  (Fig.  4.12).  For  some  type-2  traces,  there  may  be  only  one 
A  or  B  point,  or  none.  The  points  A  and  B  are  converted  to  T- 
terminations,  independently  of  each  other  and  with  probability  P[T]. 
This  procedure  introduces  some  distortion  in  the  distribution  of 
trace-2  lengths;  however,  the  resulting  bias  does  not  appear  to  be 
serious. 

The  final  result  of  applying  the  hierarchical  model  is  given  in 
Fig.  4.13,  where  a  simulation  of  Set  1  traces  through  the  midpoints  of 
Figure  4.7  and  of  Set  2  traces  using  the  conditional  model  just 
described  is  shown.  Not  only  has  each  step  in  this  model  been 
statistically  tested,  but  the  simulated  pattern  also  appears  to  be  close 
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Figure  4.12.  Galculation  of  Termination  Probability. 


to  reality  when  one  compares  Fig.  4.13  with  the  original  data  in  Fig. 
4.1b. 


4.3  Topological  Model  Applied  to  Slope  Stability  Analysis 

4.3.1  Principle 

The  hierarchical  fracture  geometry  model  was  incorporated  in 
a  topological  model  for  slope  stability  analysis.  In  its  present  state 
this  model  is  two  dimensional  like  the  hierarchical  fracture  trace 
model,  and  it  assumes  fully  persistent  fractures.  It  is  based  on  rigid 
body  assumptions.  In  this  approach  slope  stability  is  analyzed  as 
follows: 

1.  Geometric  construction  of  fracture  paths:  Depending  on  the 
actual  fracture  pattern,  either  the  homogeneous  Poisson  point  (or 
fiber)  process  model  or  the  hierarchical  fiber  process  model  can 
be  used  to  generate  the  network  of  fractures  in  a  given  two 
dimensional  slope  (the  hierarchical  fracture  model  is  at  present 
limited  to  two  dimensions). 

2.  Kinematic  analysis:  All  geometrically  possible  rock  blocks  (i.e., 
fully  persistent  blocks)  that  can  be  produced  with  the  fracture 
pattern  are  simulated  and  those  which  are  kinematically 
admissible  are  identified. 

3.  Kinetic  (or  mechanical)  analysis:  Using  an  appropriate 
mechanical  model,  the  kinetic  stability  of  the  kinematically 
admissible  blocks  is  determined, 

4.3.2  Fracture  Patterns 

The  fracture  network  is  created  using  the  hierarchical  fracture 
geometry  model.  Once  the  slope  boundaries  and  fracture  network 
are  generated,  one  needs  to  find  the  effective  fractures.  An 
"effective  fracture"  is  defined  as  a  fracture  which  can  be  a  part  (i.e.,  a 
face)  of  a  rock  block.  It  must,  therefore,  have  at  least  two 
intersection  points  with  other  fractures  or  with  the  slope  boundaries. 
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In  essence,  effective  fractures  are  inter-connected  fractures. 

Effective  fractures  are  thus  relevant  both  in  slope  stability  problems 
and  in  flow  through  rock  masses.  The  elimination  of  the  non- 
effective  fractures  which  is  done  through  sequential  iteration  is  also 
important  with  regard  to  the  computer  program  storage  design  since 
the  effective  fractures  rather  than  total  number  of  fractures  or 
intersection  points  are  used  to  define  kinematically  admissible 
fracture  paths.  As  a  consequence,  one  can  improve  the  efficiency  of 
the  path-searching  algorithm. 

A  so  called  connectivity  matrix  is  constructed  to  find  all 
fracture  paths  which  delineate  fully  persistent  rock  blocks.  Fig.  4.14 
shows  the  connectivity  matrix  and  the  searching  scheme  to  find  the 
fracture  paths.  One  essentially  starts  from  the  slope  face  and 
proceeds  from  intersection  point  to  intersection  point  to  determine 
the  fracture  paths.  These  paths  have  to  either  connect  the  face  to  the 
top  surface  or  back  to  the  face. 

4.3.3  Kinematic  Analysis 

After  constructing  the  connectivity  matrix,  all  fracture  paths 
are  checked  to  eliminate  kinematically  is  admissible  paths. 
Inadmissible  paths  are  (see  also  Figs.  4.15): 

•  The  relative  location  of  the  subsequent  intersection  point  on  a 
path  lies  in  region  IV. 

•  Parts  of  the  rock  block  are  convex  or  concave  as  shown  in  Fig. 
4.15b.  This  fact  is  recognized  in  the  search  routine  with  the 
connectivity  matrix  if  a  particular  row  in  tb  connectivity 
matrix  is  encountered  again. 

•  Face  to  face  paths  in  which  the  location  of  the  starting 
intersection  point  is  higher  than  the  end  point  (Fig.  4.15c),  even 
if  otherwise  kinematically  admissible.  This  prevents  double 
counting  of  paths. 

•  Tapered  paths  (Fig.  4.15d).  Subsequent  values  of  inclination 
angles  of  path  segments  pointing  toward  the  face  are  greater 
than  the  inclination  angle  of  the  fracture  on  which  the  block 
slides. 
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The  assumption  that  slope  failure  occurs  only  through 
movement  of  intact  blocks  and  not  by  fracturing  of  the  intact 
material  implies  that  only  sliding  along  fractures  or  separation  of 
fracture  surfaces  are  possible.  Sliding  can  kinematically  only  take 
place  along  the  fracture  with  the  smallest  dip  angle.  This  fracture 
will  be  identified  in  the  kinematic  analysis  and  so  will  be  the  ones 
which  separate. 

4.3.4  Kinetic  Analysis 

As  indicated  before,  the  kinetic  analysis  is  based  on  rigid  body 
assumptions.  Also,  only  sliding  or  separation  along  fractures  is 
allowed  but  not  toppling  and  no  fracturing  of  intact  material.  The 
kinetic  analysis  is  thus  a  limit  equilibrium  analysis  with  the  sliding 
resistance  represented  by  the  Coulomb  criterion,  i.e.  a  friction  angle 
and  a  cohesion,  and  the  assumption  of  a  so  called  "tensile  cutoff 
stress"  along  fractures  which  separate.  The  tensile  cutoff  stress 
allows  one  to  associate  a  tensile  resistance  to  fracture  separation. 

This  actually  represents  an  extension  of  the  fully  persistent  model  in 
which  such  a  resistance  would  usually  be  zero. 

4.3.5  Parametric  Studies 

Parametric  studies  were  conducted  with  the  topological  slope 
stability  model.  The  parameters  and  the  ranges  of  parameter  states 
considered  in  this  study  are  listed  in  Table  4.1.  The  effect  of  varying 
these  parameters  showed  that  the  number  of  kinematically 
admissible  rock  blocks  was  largely  controlled  by  the  number  of 
effective  fractures  which  intersect  the  slope  face,  and  that  the  most 
significant  parameters  in  these  slope  stability  analyses  were  the 
cohesion  and  friction  angle  which  controlled  the  critical  volume  and 
thus  the  probability  of  failure.  Orientation  (mean  and  dispersion) 
had  an  effect  also,  while  the  tensile  cut-off  stress  was  not  significant 
in  the  cases  studied.  Also  interesting  is  the  result  presented  in  Fig. 
4.16  which  shows  how  strongly  the  assumed  fracture  geometry 
model  affects  the  probability  of  fracturing. 


Parameter 

Fixed 

Value 

Comments 

Model 

Homogeneous  or 
Non-homogeneous 
Poisson  Model 

No.  of  Fracture 
Sets,  Nj 

1 

2 

Cohesion,  C;. 

1000 

0, 500, 
1000, 1500 

Friction 

Angle,  <j )j 

30° 

10°,  20°, 
30°,  40° 

Tensile  Cut¬ 
off  Stress,  g0j 

0 

0, 10 

50, 100 

Orientation 
Concentration 
Factor,  K 

20 

10, 20, 
30,40 

Mean  Fracture 
Orientation 
(Dip  Angle),  a 

40° 

20°,  30°, 
40°,  50° 

Mean  Trace 
Length,  m, 

6 

4, 6, 

8, 10 

Unit  Weight  of 
Rock  Mass,  y 

2200 

Fixed 

Table  4.1. 


Parameters  Considered  in  Parametric  Study 
Whenever  a  parameter  is  varied,  all  other  parametei 
assume  their  fixed  value. 
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4.4  Concluding  Comments  and  Outlook 

The  hierarchical  fracture  trace  model  produces  fracture  trace 
patterns  which  resemble  reality  very  closely.  It  represents  thus  a 
marked  improvement  over  the  previously  existing  models.  Most 
importantly,  it  is  also  capable  of  duplicating  the  sequential  fracture 
genesis  process  which  will  usually  occur  in  nature.  This  provides  a 
sound  basis  for  linking  geometric  and  mechanical  models  of  rock 
fractures  and,  very  likely  also  for  fractures  in  other  materials.  The 
topological  slope  stability  model  builds  upon  the  stochastic  fracture 
pattern  model  and  provides  a  basis  for  more  complex  stability  and 
deformation  models  of  discontinuua. 

Futuiv  wogk  in  the  stochasic  modelling  domain  will  have  to 
consider  moi^  than  two  sets  and  the  third  dimension.  The  former  is  a 
relatively  simple  extension,  the  latter  is  a  complex  problem  for  which 
a  number  of  approaches  need  to  be  investigated.  Related  to  this 
problem  is  the  consideration  of  other  fracture  characteristics  such  as 
deviation  from  a  plane  surface  and  varying  mechanical  properties. 
The  topological  model  will  have  to  be  combined  with  more 
encompassing  failure  mechanisms.  The  fracture  coalescence  model 
developed  in  this  research  and  described  in  Section  2  of  this 
Executive  Summary  (as  well  as  in  Part  I  of  the  report  series)  is  an 
obvious  choice. 
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HON- POISSON  MODEL 
HOMOGENEOUS  POISSON 


PROCESS  MODEL 


Figure  4.16.  Pf  Given  Volume  of  Rock  Block:  Variation  of  Midpoint 
Model  of  Fractures. 
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5.  Conclusions 

The  research  on  stochastic  and  centrifuge  modelling  of  jointed 
rock  resulted  in  significant  contributions  in  the  three  major  research 
areas. 

In  the  area  of  "Fracturing  of  Fractured  Rock"  we  found  that 
pre-existing  non-overlapping  fractures  coalesce  through  secondary 
cracks  and  not  through  wing  cracks.  A  predictive  model  was 
developed  and  successfully  validated.  What  has  been  discovered  and 
developed  in  this  context  does  potentially  apply  to  any  brittle 
discontinuum  and  may  thus  have  far  reaching  implications. 

In  the  area  of  centrifuge  modelling,  we  found  that  arching  in 
granular  material  under  the  realistic  geostatic  stress  conditions 
involves  large  displacements  during  which  the  load  remains  near  the 
minimum  in  contrast  to  previously  known  results  which  indicate  a 
rapid  load  increase.  An  appropriate  arching  theory  was  developed. 
This  research  also  showed  that  centrifuge  modelling  is  applicable  to 
granular  materials  with  the  usual  scaling  relationships  while 
discontinua  can  be  modelled  only  to  a  limited  extent  in  the 
centrifuge. 

In  stochastic  modelling  of  fractured  rock  we  were  able  to 
develop  a  new  hierarchical  model  which  allowed  us  to  accurately 
represent  fracture  trace  patterns  and,  most  importantly,  their 
sequential  genesis.  The  geometric  fracture  representation  was 
incorporated  in  a  topological  model  for  slope  stability  analysis. 

All  these  contributions  are  not  only  significant  with  regard  to 
the  advance  in  the  state  of  the  science  but  also  regarding  their 
practical  implications.  They  also  show  very  clearly  in  which  way  one 
should  proceed  to  produce  further  specific  and  generally  valid 
results. 


93 


6.  Acknowledgements 

The  research  results  reflect  the  active  interest  and  many 
technical  discussions  with  the  Technical  Monitors,  Dr.  S.  Wu  of  the 
U.S.  Air  Force  Office  of  Scientific  Research  and  Lt.  S.  Kuennen  of  the 
U.S.  Air  Force  Engineering  Sciences  Center.  Without  them  we  could 
not  have  achieved  what  we  did  and  this  we  would  like  to  gratefully 
acknowledge. 


94 


APPENDIX  A 


Literature  References 

ABAOUS  User's  Manual.  1990.  Hibbitt,  Karlsson  and  Sorensen,  Inc. 

Barton,  C.C.;  Larson,  E.  1985.  "Fractal  Geometry  of  Two-dimensional 
Fracture  Networks  at  Yucca  Mountain,  Southwest  Nevada,  Proc.  of 
Int.  Symp.  on  Fundamentals  of  Rock  Joints,  77-84. 

Bierbaumer,  A.  1913.  "Die  Dimensionerung  des  Tunnelmauerwerks", 
Engelmann,  Leipzig. 

Engesser,  Fr.  1882.  "Ueber  den  Erdduck  gegen  innere  Stutzwande 
(Tunnelwande)",  Deutsche  Bauzeitung.  No.  16,  pp.  91-93. 

Evans,  C.H.  1983.  "An  Examination  of  Arching  in  Granular  Soils",  S.M. 
Thesis.  Dept,  of  Civil  Eng.,  MIT. 

Iglesia,  G.;  Einstein,  H.H.;  Whitman,  R.V.;  Jessberger,  H.L.;  Guttler,  U. 
1991.  "Trapdoor  Experiments  with  Simulated  Jointed  Rock", 

Centrifuge  91.  Univ.  of  Colorado,  Boulder. 

Joseph,  P.  G.;  Einstein,  H.H.;  Whitman,  R.V.  1987.  "A  Literature  Review 
of  Geotechnical  Centrifuge  Modeling  with  Particular  Emphasis  on 
Rock  Mechanics,"  MIT  Report  to  the  Air  Force  Engineering  and 
Services  Center,  Tyndall  Air  Force  Base,  Florida. 

Joseph,  P.G.;  Einstein,  H.H.,  1987.  "Rock  Modelling  Using  the 
Centrifuge",  MIT  Report  to  the  Air  Force  Engineering  and  Services 
Center,  Tyndall  Air  Force  Base,  Florida. 

Lemaitre,  J.,  1986.  "Local  approach  of  fracture",  Engineering  Fracture 
Mechanics,  25,  5/6,  pp.  523-537. 

Reyes,  O.  M.  1988.  "Numerical  Modelling  of  Fracture  Propagation  in 
Tension  Softening  Materials",  S.M.  Thesis.  Dept,  of  Civil  Eng.,  MIT. 


Rots,  J.G.;  Nauta,  P.;  Kusters,  G.;  Blaauwendraad,  J.  1985.  "Smeared 
crack  approach  and  fracture  localization  in  concrete",  HERON.  30,  1. 


Segall,  P.  ;  Pollard,  D.D.  1983.  "Joint  Formation  in  Granitic  Rock  of  the 
Sierra  Nevada",  Geol.  Soc.  Am.  Bull.  94:563-575. 


Stone,  K.J.L.;  Guttler,  U.  1989.  "Bochum-MIT  Collaboration  on  active 
arching  experiments:  Report  from  preliminary  centrifuge  tests", 
Ruhr-University  Bochum. 

Terzaghi,  K.  1936.  "Stress  Distribution  in  Dry  and  in  Saturated  Sand 
Above  a  Yielding  Trap-Door",  Proc..  First  Int,  Conf.  on  Soil  Mech.  and 
Foundation  Eng.,  pp.  307-311. 


